Missing data, measurement error, and the future of probability

Multiple-Imputation as a way of handling Missing Data (JB)

Reading for Multiple-Imputation:

Honaker, King, and Blackwell (2011)

Bhaskaran and Smeeth (2014)

Blackwell, Honaker, and King (2017)

McElreath (2020)

https://gking.harvard.edu/category/research-interests/methods/missing-data

Missing data poses a problem. Missingness biases inference. However missingness it is not intractable problem if we assume that the mechanism giving rise to the missing data is random conditional on known features of the datset. Statistitians call this assumption “MAR: Missing at Random.”¹

Let us visualise in a subset of the the longitudional NZ dataset

First how many Id’s per wave in this dataset:

df <- nz12 %>%
  select(
    Id,
    CharityDonate,
    Emp.JobSecure,
    Male,
    Employed,
    Relid,
    Wave,
    yearS,
    KESSLER6sum,
    Age,
    yearS
  )

# always inspect your dataframe
glimpse(df)

Rows: 4,140
Columns: 10
$ Id            <fct> 15, 15, 15, 15, 15, 15, 15, 15, 15, …
$ CharityDonate <dbl> 20, 0, 5, 10, 70, 0, 170, 160, 80, 1…
$ Emp.JobSecure <dbl> 4, 6, 6, NA, 7, 5, NA, 7, NA, 7, NA,…
$ Male          <fct> Male, Male, Male, Male, Male, Male, …
$ Employed      <dbl> 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, …
$ Relid         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, …
$ Wave          <fct> 2010, 2011, 2012, 2013, 2014, 2015, …
$ yearS         <dbl> 27, 347, 834, 1200, 1608, 2037, 2336…
$ KESSLER6sum   <int> 4, 4, 4, 4, 3, 4, 5, 4, 3, 5, 1, 2, …
$ Age           <dbl> 38.82820, 39.70431, 41.03765, 42.036…

nz12 %>%
  group_by(Wave) %>%
  summarise(Unique_Id = n_distinct(Id))

# A tibble: 10 x 2
   Wave  Unique_Id
   <fct>     <int>
 1 2010        414
 2 2011        414
 3 2012        414
 4 2013        414
 5 2014        414
 6 2015        414
 7 2016        414
 8 2017        414
 9 2018        414
10 2019        414

That’s not many, but the data will be useful to explore the multiple-imputation approach

We can visualise the data, using the naniar package

library(naniar)
vis_miss(df)

We can see substantial missingness for Emp.JobSecure.

Let’s explore this:

df%>%
  select(Wave, Emp.JobSecure) %>%
  group_by(Wave)%>%
  tally(is.na(Emp.JobSecure))

# A tibble: 10 x 2
   Wave      n
   <fct> <int>
 1 2010     89
 2 2011    106
 3 2012    118
 4 2013    117
 5 2014    129
 6 2015    133
 7 2016    414
 8 2017    167
 9 2018    172
10 2019    183

Lot’s of missingness in Emp.JobSecure and the question was not included in 2016

table1::table1(~ Wave|Emp.JobSecure, data = df, overall = FALSE)

	1 (N=90)	2 (N=96)	3 (N=139)	4 (N=268)	5 (N=403)	6 (N=731)	7 (N=785)
Wave
2010	12 (13.3%)	12 (12.5%)	18 (12.9%)	31 (11.6%)	68 (16.9%)	90 (12.3%)	94 (12.0%)
2011	16 (17.8%)	14 (14.6%)	15 (10.8%)	45 (16.8%)	43 (10.7%)	84 (11.5%)	91 (11.6%)
2012	11 (12.2%)	10 (10.4%)	17 (12.2%)	41 (15.3%)	43 (10.7%)	86 (11.8%)	88 (11.2%)
2013	11 (12.2%)	11 (11.5%)	17 (12.2%)	35 (13.1%)	44 (10.9%)	90 (12.3%)	89 (11.3%)
2014	10 (11.1%)	10 (10.4%)	16 (11.5%)	28 (10.4%)	40 (9.9%)	91 (12.4%)	90 (11.5%)
2015	5 (5.6%)	10 (10.4%)	18 (12.9%)	27 (10.1%)	46 (11.4%)	78 (10.7%)	97 (12.4%)
2016	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
2017	7 (7.8%)	9 (9.4%)	12 (8.6%)	21 (7.8%)	46 (11.4%)	71 (9.7%)	81 (10.3%)
2018	8 (8.9%)	12 (12.5%)	13 (9.4%)	23 (8.6%)	35 (8.7%)	73 (10.0%)	78 (9.9%)
2019	10 (11.1%)	8 (8.3%)	13 (9.4%)	17 (6.3%)	38 (9.4%)	68 (9.3%)	77 (9.8%)

There are various methods for multiple imputation. First, let’s look at the Amelia package

library(Amelia)
# set seed
set.seed(1234)

# we need to pass a dataframe to Amelia
prep <- as.data.frame(df) # tibble won't run in amelia !!


# this is the key code
prep2 <- Amelia::amelia(
  prep,
  #dataset to impute
  m = 10,
  # number of imputations
  cs = c("Id"),
  # the cross sectional variable
  ts = c("yearS"),
  # Time series, allowing polynomials
  #ords =  none in this dataset, but use this command for ordinal data
  #logs = ,  # big numbers better to use the natural log
  sqrt = c("KESSLER6sum", "CharityDonate"),
  # skewed positive data such as K6
  noms = c("Male",  # nominal vars
           "Employed"),
  idvars = c("Wave"),
  # not imputing outcomes
  polytime = 3
) #https://stackoverflow.com/questions/56218702/missing-data-warning-r

-- Imputation 1 --

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
 61 62

-- Imputation 2 --

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 21 22 23 24 25 26 27 28

-- Imputation 3 --

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 21 22 23 24

-- Imputation 4 --

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 21 22 23

-- Imputation 5 --

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 21 22 23 24 25

-- Imputation 6 --

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
 61 62 63 64 65 66 67 68 69 70

-- Imputation 7 --

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
 41 42 43

-- Imputation 8 --

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20


-- Imputation 9 --

  1  2  3  4  5

-- Imputation 10 --

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 21 22 23 24 25 26 27 28 29 30

We can use trace plots to examine how Amelia imputed and with how much uncertainty:

Here are the imputations for Job Security (a random selection)

Amelia::tscsPlot(
  prep2,
  cs = c("15", "19", "20", "39", "549", "1078"),
  main = "Imputation of Job security",
  var = "Emp.JobSecure",
  ylim = c(0, 30)
)

Here are the imputations for Charity (another random selection). Note that there is a fair amount of uncertainty here:

Amelia::tscsPlot(
  prep2,
  cs = c("394", "1039", "1082", "340", "365", "1149", "1238" , "1253","1229"),
  main = "Impuatation of Charity",
  var = "CharityDonate",
  ylim = c(0, 10000)
)

We can center and scale our variables using the following code head(df)

prep3 <- Amelia::transform.amelia(
  prep2,
  Id = as.factor(Id),
  # redundant
  Age.10yrs = (Age / 10),
  years_s = scale(yearS, center = TRUE, scale = TRUE),
  years = yearS,
  KESSLER6sum_S = scale(KESSLER6sum, center = TRUE, scale =TRUE),
  Employed = factor(Employed),
  Relid = scale(Relid, scale = TRUE, center = TRUE),
  Male = as.factor(Male),
  Emp.JobSecure_S = scale(Emp.JobSecure, center = TRUE, scale = FALSE),
  CharityDonate = as.integer(CharityDonate)
)

# center an d scale age
out <- Amelia::transform.amelia(prep3, Age_in_Decades_C = scale(Age.10yrs,scale =FALSE, center=TRUE))

Which variables do we need to include to estimate the causal effect of job security on charity?

Write your dag!

library(ggdag)
dg <-
  dagify(
    charity ~ jobsecure + employed + age + male + relid + years,
    jobsecure ~ employed + distress + male + years + age,
    distress ~ male  + employed + years + age,
    relid ~ male + years,
    age ~ years,
    labels = c(
      "charity" = "charity",
      "jobsecure" = "job security",
      "employed" = "employed",
      "age" = "age",
      "male" = "male",
      "relid" = "religious identity",
      "years" = "years"
    ),
    exposure = "jobsecure",
    outcome = "charity"
  ) %>%
  tidy_dagitty(layout = "nicely")

ggdag::ggdag_adjustment_set(dg)

To obtain an unbiased estimate of jobsecurity on charity we must condition on employed, male, age, and years.

We can write the model using the lme4 package, which is fast. I wrote a little function, recall that we have 8 data sets

# first write out the model equation
library(lme4)
mod_eq <-  'CharityDonate  ~  Emp.JobSecure_S  + Employed + Age_in_Decades_C +  Male  +  years_s +  (1|Id)' 

# run models iterating over imputed data
loop_glmer_model <-
  function(x, y) {
    # x is the mod equation, y is the data
    m <- 10
    mod <- NULL
    for (i in 1:m) {
      mod[[i]] <- glmer(x, data = y$imputations[[i]], family = "poisson")
    }
    return(mod)
  }

m_list <- loop_glmer_model(mod_eq, out)

Here is a function for obtaining the results:

# table of effects
loop_lmer_model_tab <- function(x) {
  mp <- lapply(x, model_parameters)
  out <- parameters::pool_parameters(mp)
  return(out)
}

# create table
tab_impute <- loop_lmer_model_tab(m_list)
tab_impute

# Fixed Effects

Parameter         | Coefficient |   SE |         95% CI | Statistic |      p
----------------------------------------------------------------------------
(Intercept)       |        6.27 | 0.50 | [ 5.28,  7.26] |     12.42 | < .001
Emp.JobSecure_S   |        0.03 | 0.03 | [-0.03,  0.08] |      0.95 | 0.342 
Employed [1]      |        0.16 | 0.03 | [ 0.09,  0.23] |      4.56 | < .001
Age_in_Decades_C  |       -4.41 | 1.27 | [-6.90, -1.92] |     -3.47 | < .001
Male [Not_Male]   |       -1.33 | 0.69 | [-2.69,  0.02] |     -1.92 | 0.054 
years_s           |        1.26 | 0.37 | [ 0.53,  1.98] |      3.41 | < .001
SD (Intercept)    |        6.31 |      |                |           |       
SD (Observations) |        1.00 |      |                |           |       

# create graph
plot_impute <- plot(tab_impute)
plot_impute

We can plot the effects using a coefficient plot

library(ggeffects)
library(gghighlight) # not used here, useful for interactions 
graph_predictions_imputed <- function(x, y) {
  # x = model objects
  m <- 10
  out <- NULL
  for (i in 1:m) {
    out[[i]] <-
      ggeffects::ggpredict(x[[i]], terms = c("Emp.JobSecure_S"))
  }
  plots <- NULL
  for (i in 1:m) {
    plots[[i]] <-
      plot(out[[i]], facets = T) # + scale_y_continuous(limits=c(6.35,6.85) )
  }
  plots[[10]] +
    gghighlight::gghighlight() +
    ggtitle(y)
}

# graph
graph_predictions_imputed(m_list,"Effect of Jobsecurity on Charity (not reliable")

If you want a LaTeX table, you can use this code:

library(huxtable)

huxtable::as_hux( your_model_here ) %>%
  select("Parameter", "Coefficient", "CI_low", "CI_high", "p") %>%
  set_number_format(3) %>%
  set_left_padding(20) %>%
  set_bold(1, everywhere) %>%
  quick_latex()

Compare imputation results with row-wise deleted results

When you run a regression with missing data, R automateically deletes the missing cases.

Let’s look at the results from the row-wise deleted data:

# prepare data as we did for the imputated dataset
df2 <- df %>%
  dplyr::mutate(
    Age.10yrs = (Age / 10),
    Age_in_Decades_C = scale(Age.10yrs, scale = FALSE, center = TRUE),
    years_s = scale(yearS, center = TRUE, scale = TRUE),
    years = yearS,
    KESSLER6sum_S = scale(KESSLER6sum, center = TRUE, scale = TRUE),
    Employed = factor(Employed),
    Relid = scale(Relid, scale = TRUE, center = TRUE),
    Male = as.factor(Male),
    Emp.JobSecure_S = scale(Emp.JobSecure, center = TRUE, scale = FALSE)
  )

# run model
m_no_impute <- glmer(mod_eq, data = df2, family = "poisson")

# create table
tab_no <-
  parameters::model_parameters(m_no_impute, effects = c("all"))
tab_no

# Fixed Effects

Parameter        |  Log-Mean |       SE |         95% CI |      z |      p
--------------------------------------------------------------------------
(Intercept)      |      5.26 |     0.08 | [ 5.09,  5.42] |  62.00 | < .001
Emp.JobSecure_S  | -5.65e-03 | 5.87e-04 | [-0.01,  0.00] |  -9.61 | < .001
Employed [1]     |      0.19 | 6.46e-03 | [ 0.18,  0.20] |  29.18 | < .001
Age_in_Decades_C |     -0.89 |     0.04 | [-0.96, -0.81] | -21.86 | < .001
Male [Not_Male]  |     -0.64 |     0.11 | [-0.86, -0.43] |  -5.79 | < .001
years_s          |      0.18 |     0.01 | [ 0.16,  0.21] |  15.41 | < .001

# Random Effects

Parameter          | Coefficient
--------------------------------
SD (Intercept: Id) |        1.00
SD (Residual)      |        1.00

# create graph
plot_no <- plot(tab_no)
plot_no

When we compare the graphs, we see that the multiply imputed datasets shrink estimates towards zero.

Multiple imputation is sometimes avoided because people don’t like to “invent” data. However, creating multiply imputed datasets and integrating over their uncertainty during model tends to increase uncertainty in a model. That’s generally a good thing when we want to predict features of the population.

library(patchwork)

plot_impute / plot_no +
  plot_annotation(title = "Comparison of regressions using (a) multiple-imputed  and (b) row-wise deleted datasets",
                  tag_levels = 'a')

However it would be a mistake to think that multiple imputation is sufficient. Note that the case-wise deleted data is confident that men give less to charity than do women and other genders. However the model is tending to inflate estimates for men, perhaps because men tend to disproportionately leave this question unanswered. Let’s check this intuition:

df %>%
  dplyr::mutate(CharityNa = is.na(CharityDonate)) %>%
  count(Male, CharityNa) %>%
  group_by(Male) %>%
  mutate(freq = n / sum(n))

# A tibble: 4 x 4
# Groups:   Male [2]
  Male     CharityNa     n   freq
  <fct>    <lgl>     <int>  <dbl>
1 Male     FALSE      1506 0.972 
2 Male     TRUE         44 0.0284
3 Not_Male FALSE      2462 0.951 
4 Not_Male TRUE        128 0.0494

No, the intution was wrong. However, it is important to realise that missingness might not be completely at random conditional on variables in the model. In this case, you dataset cannot avoid its biases by multiple imputation.

Recall that it is generally useful to include more variables for accurate prediction. For this reason, when multiply imputing outcomes it is generally useful to include more variables.

Appendix 1: Imputation using BRMS

Multiple imputation in BRMS (warning: these models take 10x longer)

# prepare data
# you must create a list of the imputed datasets, like so:

brmdat <- list(
  out$imputations$imp1,
  
  out$imputations$imp2,
  
  out$imputations$imp3,
  
  out$imputations$imp4,
  
  out$imputations$imp5,
  
  out$imputations$imp6,
  
  out$imputations$imp7,
  
  out$imputations$imp8,
  
  out$imputations$imp9,
  
  out$imputations$imp10
)

We write the model as follows:

# model equation
# note the `bf` command

bf_mod_eq <- bf(mod_eq)

# we use the `brm_multiple` syntax and feed in the brmdat

fit_imp1 <-
  brm_multiple(
    bf_mod_eq,
    data = brmdat,
    family = "poisson",
    file = here::here("models", "bayes-imp-1")
  )

# table
summary(fit_imp1)

 Family: poisson 
  Links: mu = log 
Formula: CharityDonate ~ Emp.JobSecure_S + Employed + Age_in_Decades_C + Male + years_s + (1 | Id) 
   Data: brmdat (Number of observations: 4140) 
Samples: 40 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 40000

Group-Level Effects: 
~Id (Number of levels: 414) 
              Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)     5.40      2.89     1.08     9.93 7.38
              Bulk_ESS Tail_ESS
sd(Intercept)       41       47

Population-Level Effects: 
                 Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept            7.33      1.83     4.51    11.56 8.47
Emp.JobSecure_S      0.03      0.02    -0.01     0.06 5.11
Employed1            0.16      0.03     0.11     0.20 4.43
Age_in_Decades_C    -3.14      2.04    -5.65     0.13 5.31
MaleNot_Male        -0.55      0.87    -2.81     0.56 7.08
years_s              0.89      0.59    -0.06     1.62 5.36
                 Bulk_ESS Tail_ESS
Intercept              41       40
Emp.JobSecure_S        42       86
Employed1              43       68
Age_in_Decades_C       42       49
MaleNot_Male           41       48
years_s                42       53

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit_imp1)  # VERY POOR MIXING

Mixing is a problem, indeed I don’t think I have ever run a model that has mixed worse than this one.

We see the trouble in the coefficient plots:

#plot posteriors
plot_bayes_1 <- brms::mcmc_plot(fit_imp1,
                                type = "areas",
                                prob = .89)

plot_bayes_1 +
  plot_impute +
  plot_annotation(title = "Comparison of regressions using  (a) Bayesian Imputationand and (b) Lmer models",
                  tag_levels = 'a')

We can model missing-ness in BRMS in one step (for continuous missing variables). Let’s try that next:

# Emp.JobSecure_S  + Employed + Age_in_Decades_C +  Male  +  years_s

#Note that BRMS can only impute continuous data. That's not a problem. Your factors are converted to integars anyway.
df3 <- df2%>%
  dplyr::mutate(Employed = as.numeric(as.character(Employed)),
         Male = as.numeric((Male))-1) # male as zero or 1

glimpse(df3)

Rows: 4,140
Columns: 16
$ Id               <fct> 15, 15, 15, 15, 15, 15, 15, 15, 1…
$ CharityDonate    <dbl> 20, 0, 5, 10, 70, 0, 170, 160, 80…
$ Emp.JobSecure    <dbl> 4, 6, 6, NA, 7, 5, NA, 7, NA, 7, …
$ Male             <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, …
$ Employed         <dbl> 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, …
$ Relid            <dbl[,1]> <matrix[20 x 1]>
$ Wave             <fct> 2010, 2011, 2012, 2013, 2014,…
$ yearS            <dbl> 27, 347, 834, 1200, 1608, 2037, 2…
$ KESSLER6sum      <int> 4, 4, 4, 4, 3, 4, 5, 4, 3, 5, 1, …
$ Age              <dbl> 38.82820, 39.70431, 41.03765, 42.…
$ Age.10yrs        <dbl> 3.882820, 3.970431, 4.103765, 4.2…
$ Age_in_Decades_C <dbl[,1]> <matrix[20 x 1]>
$ years_s          <dbl[,1]> <matrix[20 x 1]>
$ years            <dbl> 27, 347, 834, 1200, 1608, 2037, 2…
$ KESSLER6sum_S    <dbl[,1]> <matrix[20 x 1]>
$ Emp.JobSecure_S  <dbl[,1]> <matrix[20 x 1]>

df3%>%
  dplyr::count(CharityDonate ==0)

# A tibble: 3 x 2
  `CharityDonate == 0`     n
  <lgl>                <int>
1 FALSE                 3453
2 TRUE                   515
3 NA                     172

## Write the model, note the `mi`s: each `mi` needs to appear as the outcome of a model. 
## note that for simplicity, these models claim that mi is random conditional on year and individual. We might obtain better predictions of missingness by including more information. 

# Note that brms can only take identity link functions when handling missing data. The "lognormal" is better than the normal because the variances are estimated using a log link. 

bform <- 
  bf(CharityDonate + 1 | mi() ~ mi(Emp.JobSecure_S) + mi(Employed) + mi(Age_in_Decades_C)  + mi(Male)  + (1 + years_s | Id), family = "lognormal") +
  bf(Emp.JobSecure_S | mi() ~ mi(Employed) + mi(Age_in_Decades_C)  + mi(Male)  + (1 + years_s| Id)) +
  bf(Employed | mi() ~ (1 + years_s| Id)) +
  bf(Age_in_Decades_C | mi() ~  (1 + years_s| Id))+
  bf(Male | mi() ~ (1 + years_s | Id)) +
  set_rescor(FALSE)

## fit the model 
fit_imp2 <- brm(bform, 
                data = df3, 
                file = here::here("models", "bayes-imp-2"))
summary(fit_imp2)

 Family: MV(lognormal, gaussian, gaussian, gaussian, gaussian) 
  Links: mu = identity; sigma = identity
         mu = identity; sigma = identity
         mu = identity; sigma = identity
         mu = identity; sigma = identity
         mu = identity; sigma = identity 
Formula: CharityDonate + 1 | mi() ~ mi(Emp.JobSecure_S) + mi(Employed) + mi(Age_in_Decades_C) + mi(Male) + (1 + years_s | Id) 
         Emp.JobSecure_S | mi() ~ mi(Employed) + mi(Age_in_Decades_C) + mi(Male) + (1 + years_s | Id) 
         Employed | mi() ~ (1 + years_s | Id) 
         Age_in_Decades_C | mi() ~ (1 + years_s | Id) 
         Male | mi() ~ (1 + years_s | Id) 
   Data: df3 (Number of observations: 4140) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~Id (Number of levels: 414) 
                                                     Estimate
sd(CharityDonate1_Intercept)                             1.25
sd(CharityDonate1_years_s)                               0.16
sd(EmpJobSecureS_Intercept)                              0.35
sd(EmpJobSecureS_years_s)                                0.17
sd(Employed_Intercept)                                   0.13
sd(Employed_years_s)                                     0.04
sd(AgeinDecadesC_Intercept)                              0.55
sd(AgeinDecadesC_years_s)                                0.07
sd(Male_Intercept)                                       0.36
sd(Male_years_s)                                         0.00
cor(CharityDonate1_Intercept,CharityDonate1_years_s)    -0.36
cor(EmpJobSecureS_Intercept,EmpJobSecureS_years_s)       0.25
cor(Employed_Intercept,Employed_years_s)                 0.15
cor(AgeinDecadesC_Intercept,AgeinDecadesC_years_s)      -0.30
cor(Male_Intercept,Male_years_s)                         0.26
                                                     Est.Error
sd(CharityDonate1_Intercept)                              0.70
sd(CharityDonate1_years_s)                                0.19
sd(EmpJobSecureS_Intercept)                               0.46
sd(EmpJobSecureS_years_s)                                 0.18
sd(Employed_Intercept)                                    0.17
sd(Employed_years_s)                                      0.06
sd(AgeinDecadesC_Intercept)                               0.32
sd(AgeinDecadesC_years_s)                                 0.09
sd(Male_Intercept)                                        0.10
sd(Male_years_s)                                          0.00
cor(CharityDonate1_Intercept,CharityDonate1_years_s)      0.56
cor(EmpJobSecureS_Intercept,EmpJobSecureS_years_s)        0.55
cor(Employed_Intercept,Employed_years_s)                  0.31
cor(AgeinDecadesC_Intercept,AgeinDecadesC_years_s)        0.26
cor(Male_Intercept,Male_years_s)                          0.31
                                                     l-95% CI
sd(CharityDonate1_Intercept)                             0.07
sd(CharityDonate1_years_s)                               0.00
sd(EmpJobSecureS_Intercept)                              0.00
sd(EmpJobSecureS_years_s)                                0.01
sd(Employed_Intercept)                                   0.00
sd(Employed_years_s)                                     0.00
sd(AgeinDecadesC_Intercept)                              0.04
sd(AgeinDecadesC_years_s)                                0.00
sd(Male_Intercept)                                       0.20
sd(Male_years_s)                                         0.00
cor(CharityDonate1_Intercept,CharityDonate1_years_s)    -0.97
cor(EmpJobSecureS_Intercept,EmpJobSecureS_years_s)      -0.46
cor(Employed_Intercept,Employed_years_s)                -0.35
cor(AgeinDecadesC_Intercept,AgeinDecadesC_years_s)      -0.72
cor(Male_Intercept,Male_years_s)                        -0.22
                                                     u-95% CI
sd(CharityDonate1_Intercept)                             1.89
sd(CharityDonate1_years_s)                               0.47
sd(EmpJobSecureS_Intercept)                              1.15
sd(EmpJobSecureS_years_s)                                0.45
sd(Employed_Intercept)                                   0.41
sd(Employed_years_s)                                     0.13
sd(AgeinDecadesC_Intercept)                              0.89
sd(AgeinDecadesC_years_s)                                0.23
sd(Male_Intercept)                                       0.45
sd(Male_years_s)                                         0.00
cor(CharityDonate1_Intercept,CharityDonate1_years_s)     0.33
cor(EmpJobSecureS_Intercept,EmpJobSecureS_years_s)       0.97
cor(Employed_Intercept,Employed_years_s)                 0.47
cor(AgeinDecadesC_Intercept,AgeinDecadesC_years_s)       0.02
cor(Male_Intercept,Male_years_s)                         0.65
                                                     Rhat
sd(CharityDonate1_Intercept)                         3.31
sd(CharityDonate1_years_s)                           3.47
sd(EmpJobSecureS_Intercept)                          3.78
sd(EmpJobSecureS_years_s)                            3.10
sd(Employed_Intercept)                               3.01
sd(Employed_years_s)                                 3.48
sd(AgeinDecadesC_Intercept)                          3.31
sd(AgeinDecadesC_years_s)                            2.94
sd(Male_Intercept)                                   3.57
sd(Male_years_s)                                     4.27
cor(CharityDonate1_Intercept,CharityDonate1_years_s) 4.06
cor(EmpJobSecureS_Intercept,EmpJobSecureS_years_s)   3.44
cor(Employed_Intercept,Employed_years_s)             3.12
cor(AgeinDecadesC_Intercept,AgeinDecadesC_years_s)   3.27
cor(Male_Intercept,Male_years_s)                     3.42
                                                     Bulk_ESS
sd(CharityDonate1_Intercept)                                4
sd(CharityDonate1_years_s)                                  4
sd(EmpJobSecureS_Intercept)                                 4
sd(EmpJobSecureS_years_s)                                   5
sd(Employed_Intercept)                                      5
sd(Employed_years_s)                                        4
sd(AgeinDecadesC_Intercept)                                 4
sd(AgeinDecadesC_years_s)                                   5
sd(Male_Intercept)                                          4
sd(Male_years_s)                                            4
cor(CharityDonate1_Intercept,CharityDonate1_years_s)        4
cor(EmpJobSecureS_Intercept,EmpJobSecureS_years_s)          4
cor(Employed_Intercept,Employed_years_s)                    5
cor(AgeinDecadesC_Intercept,AgeinDecadesC_years_s)          4
cor(Male_Intercept,Male_years_s)                            4
                                                     Tail_ESS
sd(CharityDonate1_Intercept)                               11
sd(CharityDonate1_years_s)                                 16
sd(EmpJobSecureS_Intercept)                                11
sd(EmpJobSecureS_years_s)                                  18
sd(Employed_Intercept)                                     20
sd(Employed_years_s)                                       19
sd(AgeinDecadesC_Intercept)                                15
sd(AgeinDecadesC_years_s)                                  32
sd(Male_Intercept)                                         12
sd(Male_years_s)                                           11
cor(CharityDonate1_Intercept,CharityDonate1_years_s)       14
cor(EmpJobSecureS_Intercept,EmpJobSecureS_years_s)         13
cor(Employed_Intercept,Employed_years_s)                   15
cor(AgeinDecadesC_Intercept,AgeinDecadesC_years_s)         23
cor(Male_Intercept,Male_years_s)                           15

Population-Level Effects: 
                                  Estimate Est.Error
CharityDonate1_Intercept              4.46      0.03
EmpJobSecureS_Intercept               0.11      0.12
Employed_Intercept                    0.70      0.01
AgeinDecadesC_Intercept              -0.02      0.05
Male_Intercept                        0.63      0.01
CharityDonate1_miEmp.JobSecure_S      0.01      0.02
CharityDonate1_miEmployed             0.35      0.10
CharityDonate1_miAge_in_Decades_C     0.15      0.06
CharityDonate1_miMale                -0.15      0.13
EmpJobSecureS_miEmployed              0.08      0.07
EmpJobSecureS_miAge_in_Decades_C      0.09      0.05
EmpJobSecureS_miMale                 -0.25      0.13
                                  l-95% CI u-95% CI Rhat
CharityDonate1_Intercept              4.43     4.50 3.22
EmpJobSecureS_Intercept              -0.01     0.29 3.30
Employed_Intercept                    0.69     0.72 3.49
AgeinDecadesC_Intercept              -0.08     0.04 3.30
Male_Intercept                        0.62     0.65 3.29
CharityDonate1_miEmp.JobSecure_S     -0.02     0.04 3.90
CharityDonate1_miEmployed             0.21     0.48 3.34
CharityDonate1_miAge_in_Decades_C     0.06     0.21 3.01
CharityDonate1_miMale                -0.30     0.05 3.43
EmpJobSecureS_miEmployed             -0.03     0.17 3.98
EmpJobSecureS_miAge_in_Decades_C      0.04     0.17 3.54
EmpJobSecureS_miMale                 -0.45    -0.10 3.34
                                  Bulk_ESS Tail_ESS
CharityDonate1_Intercept                 4       15
EmpJobSecureS_Intercept                  4       19
Employed_Intercept                       4       16
AgeinDecadesC_Intercept                  4       12
Male_Intercept                           4       14
CharityDonate1_miEmp.JobSecure_S         4       12
CharityDonate1_miEmployed                4       11
CharityDonate1_miAge_in_Decades_C        5       22
CharityDonate1_miMale                    4       15
EmpJobSecureS_miEmployed                 4       12
EmpJobSecureS_miAge_in_Decades_C         4       15
EmpJobSecureS_miMale                     4       13

Family Specific Parameters: 
                     Estimate Est.Error l-95% CI u-95% CI
sigma_CharityDonate1     1.69      0.45     1.34     2.47
sigma_EmpJobSecureS      1.43      0.20     1.09     1.63
sigma_Employed           0.39      0.10     0.21     0.46
sigma_AgeinDecadesC      0.74      0.60     0.01     1.39
sigma_Male               0.00      0.00     0.00     0.00
                     Rhat Bulk_ESS Tail_ESS
sigma_CharityDonate1 3.45        4       16
sigma_EmpJobSecureS  3.12        5       11
sigma_Employed       3.64        4       11
sigma_AgeinDecadesC  3.09        5       20
sigma_Male           3.22        4       11

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Again, the model is a all over the place:

plot(fit_imp2)

Mixing was a problem. Minimally we’d need to include more variables. However there might be deeper problems with the data, or our approach. Analysis is iterative. We need to return to the data and scrutinise it more carefully.

plot_bayes_2 <- brms::mcmc_plot(fit_imp2,
                                type = "areas",
                                prob = .89)

plot_bayes_2

# again we run into trouble. 
plot(fit_imp2)

Appendix 2. Probability

We have been working with probability throughout this course.

Suppose there is a test that is 99% accurate at detecting COVID if you have it.

Very rarely it throws up a false positive,say one in a thousand.

You just tested positve. What is the probability that you have COVID? Our intuition is that we probably have COVID. However, let’s assume COVID is rare. Currently in NZ, there are about 50 cases, so 1 in 100,000. The background rate matters.

Bayes rule says

\[ Pr(COVID|Positive) = \frac{Pr(Positive|COVID)\times Pr (COVID}{Pr(Positive)} \]

We plug in the numbers:

Pr_Positive_COVID <- 0.99
Pr_Positive_Healthy <- 0.01
Pr_COVID <- 0.00001

# Calculate the background probability of testing positive

Pr_Positive <- Pr_Positive_COVID * Pr_COVID + Pr_Positive_Healthy * ( 1 - Pr_COVID )

# Now calculated your probability of testing positive

Pr_COVID_Positive <- Pr_Positive_COVID * Pr_COVID / (Pr_Positive )
Pr_COVID_Positive # 1 in 1000

[1] 0.0009890307

---

1. For a wonderful explanation, see: here↩︎